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ABSTRACT 


Blind predictions were conducted to validate a discrete crack methodology based 
on the Floating Node Method to simulate matrix-crack/delamination interaction. The 
main novel aspects of the approach are: (1) the implementation of the floating node 
method via an “extended interface element” to represent delaminations, matrix-cracks 
and their interaction, (2) application of directional cohesive elements to infer overall 
delamination direction, and (3) use of delamination direction and stress state at the 
delamination front to determine migration onset. 

Overall, good agreement was obtained between simulations and experiments. 
However, the validation exercise revealed the strong dependence of the simulation of 
matrix-crack/delamination interaction on the strength data (in this case transverse 
interlaminar strength, Yr) used within the cohesive zone approach applied in this 
work. This strength value, Yr, is itself dependent on the test geometry from which the 
strength measurement is taken. Thus, choosing an appropriate strength value becomes 
an ad-hoc step. As a consequence, further work is needed to adequately characterize 
and assess the accuracy and adequacy of cohesive zone approaches to model small 
crack growth and crack onset. Additionally, often when simulating damage 
progression with cohesive zone elements, the strength is lowered while keeping the 
fracture toughness constant to enable the use of coarser meshes. Results from the 
present study suggest that this approach is not recommended for any problem 
involving crack initiation, small crack growth or multiple crack interaction. 
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INTRODUCTION 


As part of the NASA Advanced Composites Project (ACP), the clamped tapered 
beam specimen, illustrated in Figure 1, has been designed to provide validation data 
for progressive damage analysis models [1]. The present work reports on blind 
predictions performed with a 3D discrete crack approach based on the Floating Node 
Method (FNM) [2], of the quasi-static simulation of the clamped tapered beam. 

The Clamped Tapered Beam (CTB) specimen is based on the delamination 
migration test proposed in [3, 4]. The tapered geometry was devised to localize the 
first damage occurrence in the tapered region, without prescribing an initial crack. The 
specimen is clamped on both ends and loaded via a pin, connected to a loading rod, 
located at a distance La or Le from the left-hand-side clamp, as shown in Figure 1. The 
boundary and loading conditions were chosen to favor delamination growth and 
subsequent migration after a certain amount of delamination growth. The typical 
sequence of events consists of the onset of a matrix-crack at the tapered region which 
triggers delamination, followed by subsequent delamination migration to a different 
interface via a dominant matrix-crack. Varying the load-offset application point from 
La to Le affects both the peak load and the migration location. The specimen was 
designed to ensure that the fatigue crack growth would reside in the Paris Law crack 
growth regime of IM7/8552 Carbon Fiber Reinforced Plastic (CFRP) for the duration 
of a complete set of damage events (matrix-crack initiation, delamination, migration, 
etc.), thus avoiding premature damage arrestment or unstable failure. Further details of 
the specimen and design considerations are given in [1]. 

In the present work an approach based on the FNM is assessed. The FNM 
represents discrete crack networks by generating sub-elements within cracked finite 
elements through the use of floating nodes [2]. No a priori location for damage 
initiation is assumed. As proposed, the methodology does not limit the number of 
possible matrix-cracks, delaminations, or the crack spacing, thereby enabling the 
model to accommodate the initiation and propagation of complex crack networks at 
the meso-scale without the need for re-meshing. For models under quasi-static 
loading, the FNM approach is combined with Directional Cohesive Zone Elements 
(DCZE) to model onset and growth of delamination and matrix-cracks. 
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Figure 1. Clamped tapered beam specimen and boundary conditions. 


METHODOLGY 


In this section the methodology used in the blind-predictions exercise will be 
outlined, with emphasis on the extended interface floating node element developed, 
the directional cohesive zone elements, and the migration onset criterion. 


Floating Node Method extended interface element 


An extended interface element has been used to capture the damage morphology 
in the central region of the CTB specimen. It consists of a 3D 48-node element, 
composed of 16 nodes (named real) which have pre-assigned positions, and 32 
floating nodes (Figure 2a). The real nodes are used to define two sub-elements: one 
above and another below a pre-assigned cohesive element located at the interface 
between the two sub-elements - hence the name “extended interface” element. The 32 
floating nodes are the minimum subset of floating nodes that enable the two sub- 
elements to split independently. The floating nodes have their topological relationships 
to the element assigned in the same manner as regular nodes, and are included in the 
element connectivity list. However, if the element is not split, their degrees of freedom 
are not activated, and they can be removed from the system of equations. If the sub- 
elements need to be split to accommodate matrix-cracks, the floating nodes are used as 
required. Since the connectivity of the floating nodes is defined ‘a priori’, when 
activated, they automatically enforce crack path continuity. For example, an edge will 
be split by a unique pair of floating nodes shared by all the elements sharing that edge, 
as defined by their connectivity list. The element, as proposed, can be used to 
represent matrix-cracks with any orientation in the xy plane (Figure 2). As matrix- 
cracks initiate, Figure 2b, the interface element may also split using floating nodes, 
and cohesion elements are assigned as required to model the subsequent opening of 
the newly generated splits. The correct representation of the kinematics has been 
demonstrated to be an important feature to correctly capture matrix- 
crack/delamination interaction [2, 5]. 
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(a) Pristine extended interface element with sub-elements (b) element representing a single 
above and below a Directional Cohesive Zone Element matrix-crack 
(DCZE). 
Figure 2. Extended interface element. Distances between between sub-elements and partitions are 
exaggerated for illustration purposes. 


Matrix-cracks and through-thickness straight crack approximations 


As depicted in Figure 2b, the matrix-cracks are assumed to be straight through the 
ply thickness. This approximation is used to facilitate the integration of partitioned 
elements. However, matrix-cracks will only be straight through-thickness, as depicted 
in Figure 2b, if forming under transverse tensile loads; otherwise their angle may vary, 
see Figure 3. In order to attenuate the effect of this approximation, the present study 
uses two corrections. The first consists in determining the cohesive tractions in a 
rotated coordinate system, following [6, 7]. In addition, since the matrix-cracks are 
free to seek a mode I path as they propagate through-thickness, the shear component 
determined along direction | in Figure 3 is ignored, if not zero, and not considered to 
contribute to the mode-mixity. 


Figure 3. Matrix-crack through thickness approximation. 


Matrix-cracks are considered to onset via the Benzeggagh-Kenane (BK) criterion, 
written in traction space [8], when the effective traction: 


(Te)” = ((03))? + (sh )7B" (1) 
exceeds the strength, given by: 


(Yu)? = (Yr)? + (% 7B" (2) 


where T.;, = / (11)? + (Tz)? and (x) is defined as (x) = a(x + |x|), and Y; and Y; 


designate the transverse and shear strength, respectively; 7 is a fitting parameter; and 
B is the mode-mixity which is assumed, for matrix-cracks, to be given, by: 


(03) 
= 3 
t= eye Gay 2 


The tractions are determined in the plane perpendicular to the maximum principal 
stress obtained in the ‘mn’ plane transverse to the fibers (hence T; = 0), Figure 3. In 
addition, matrix-cracks can also be onset via the migration criterion that will be 
outlined later in this section. As mentioned previously, once a matrix-crack onsets and 
a given element is split, a cohesive element is automatically inserted and, if sufficient 


traction is applied, it will subsequently open according to its traction-separation law. 
The mixed-mode traction separation law implemented is based on the formulation 
proposed in [8]. Using Eqs. 1 and 2 to determine matrix-crack onset ensures 
consistency throughout the crack onset, insertion and subsequent opening process, 
since Eqs. 1 and 2 are also used in [8] to determine the peak of the bi-linear traction- 
separation law. 


Directional cohesive elements - active ply determination 


Directional cohesive elements are used to estimate the delamination growth 
direction. This information is in turn used to determine the ply towards which 
delamination tends to grow, here named as the “activated ply” [9]. We assume the 
crack growth direction vector d is opposite to the gradient of the displacement jump 
6, determined within the element: 


d = -V6, (4) 


where: 


bp = [52 + 63 + (53 )? (5) 


and 6; are the shear and normal displacement jump components. The gradient V6 is 
computed at the centroid, x,, of the element using the element shape functions: 
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In the present work, delamination growth direction is used in conjunction with shear 
traction vector T: 
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to determine the activated ply. This procedure is illustrated in Figure 4 and it assumes 
that microcracks form ahead of the delamination, with an orientation governed by the 
principle tensile stress. Therefore, knowing the crack direction, d, and the shear 
traction direction, t, the activated ply can be determined by: 


t-d> 0.0 : lower ply is activated 
ie : ee (8 
t-d< 0.0 :upper ply is activated 
as graphically illustrated in Figure 4. If t-d is positive, the ply below the interface 
element is activated, otherwise the ply above is considered to be activated. If t- d = 
0, both plies are considered to be activated. The criterion above assumes a stress 
normal to the interface, o > 0. If o < 0, the orientation of the matrix-cracks will 


rotate, and the criterion will be inversed, i.e., if t.d is positive, the ply above the 
interface is activated, and if t. d is negative the ply below is activated. 

In the present work, the information regarding an activated ply is used to 
determine propensity for migration. It is envisaged that in future studies this 
information can also be used to adjust for different mixed-mode fracture toughness 
values as a function of activated ply, e.g., in tape/fabric interfaces [10]. 

The cohesive zone elements’ constitutive response follows [8]. One aspect worth 
noting of such formulation is the requirement that: 


(9) 


where G;, and G,;, are the Mode I and II fracture toughness, respectively, which 
ensures that large crack growth under mode-mixed conditions is accurately predicted. 
However, as a consequence, Eq. 9 may limit the ability of the approach to accurately 
capture damage initiation and small crack growth, since the actual material properties 
may not relate as given in Eq. 9. 
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Figure 4. Determination of activated ply as a function of crack and shear traction vector direction. 


Migration onset 


Migration onset in the CTB specimen is defined as the onset of a dominant matrix- 
crack that kinks out of the plane of its originating delamination, towards a different ply 
interface. In the current analysis, migration onset is modeled by splitting elements and 
inserting cohesive elements along the predicted crack surface (or path). The criterion 
for migration onset consists of comparing microcrack growth direction to the fiber 
direction of an active ply, f,. Since it is known that the in-plane orientation of the 
microcracks is perpendicular to the shear traction vector (in-plane), the microcracks 
are assumed to have a propensity for meandering through the active ply if: 


It: f,| < COS Vor (10) 


where y,, is used to define a range of relative orientations between tT and f, that may 
lead to migration. This criterion is illustrated in Figure 5, for the case depicted in 
Figure 4. Figure 5a illustrates a case in which the criterion proposed in Eq. 10 is not 
met. Indeed, for migration onset to occur in this case, the microcracks would need to 
either break fibers or rotate by 90° as they grow into the activated ply. Figure 5b 
represents a case where Eq. 10 is met. In this case, the relative orientation of the 
matrix-cracks to the activated ply fiber direction is such that the microcracks are 
assumed to be able to meander through the activated ply. Although the value for y,,. = 
30° used in this study is based on the conditions for containment observed in [11], Yer 
is considered to be a numerical parameter that serves to prevent spurious insertion of 
matrix-cracks, while guaranteeing the delamination and migration onset are treated as 
part of the same process, as observed experimentally. The condition above is assessed 
as the traction at the cohesive element located at the interface equals the strength. If 
the criterion is met, cohesive elements are inserted through-thickness, splitting the 
sub-element. The opening of these cohesive elements, as well as the subsequent 
opening of the also-split interface cohesive element is governed by their traction- 
separation law. Hence, migration will only be completed if it is energetically 
favorable. During that assessment, the in-plane failure criterion is deactivated in the 
element being assessed. Once migration onsets, it is assumed to grow through the 
thickness with an angle of 8,,, = 45° relative to the delamination growth direction. 
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(b) In-plane view of activated ply of Figure 4b. 
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(a) In-plane view of activated ply of Figure 4a. 


Figure. 5 Illustration of migration onset criterion, Eq. 10. The criterion is met in (b) and not met in (a). In 
(a) the microcracks would need to either break fibers or rotate by 90° to grow into the activated ply, 
while in (b) the microcracks are able to meander through the activated ply with less changes to their 
path. 


NUMERICAL MODEL AND MATERIAL PROPERTIES 


The boundary conditions for the model were chosen based on the best fit to Digital 
Image Correlation (DIC) data. Further details on this procedure are given in [1]. The 
boundary conditions consist of pressure applied at the clamp region in a first step, 
followed by displacement applied at a line of nodes along the hinge as illustrated in 
Figure 1. The material properties used are summarized in Tables I and IL. Two values 
of Yr are reported in the table. These two values define lower and upper bounds of 
material strength values as measured by 90’ tensile or 3-point bend tests, and will be 
used to assess the sensitivity of the solution on Yr. 


TABLE I. ELASTIC PROPERTIES [12]. 


Ey4 E22 E33 Gy. G43 G23 Vi20 V13——«V23 
[GPa] [GPa] [GPa] [GPa] [GPa] [GPa] 
jbo yee: 8.96 8.96 5.08 5.08 2.99 0.32 0.32 0.5 
Yr Ys Gic Gic 
[MPa] [MPa] [N.mm!] [N.mm!] 
{64 [13], 127 [14]} {112,223} 0.24[15] 0.739 [16] 2.1 [17] 


Following [8], the values for Y; are assumed to be given by Eq. 9. As mentioned 
previously, while Eq. 9 guarantees the cohesive zone approach implemented is 
capable of modeling large crack propagation under mixed-mode conditions, it may 
limit its accuracy in the modeling of crack initiation and small crack growth, 
particularly for shear dominated cracks, since the actual material properties may not 
relate as given in Eq. 9. Intra-laminar and inter-laminar strength and toughness 
properties are assumed to be the same. 

The finite element mesh used for the blind-predictions is illustrated in Figure 6. 
The floating node extended interface element was used to model damage onset and 
propagation in the center region of the laminate, (Figure 6). The remainder of the 
specimen was modeled with a combination of solid (C3D8I) and shell (SC8R) 
Abaqus® elements. A mesh refinement study was performed for both Yt = 64 MPa 
and Yr = 127 MPa, focusing on mesh convergence of the peak load value. This study 
showed that the peak load predictions obtained by halving the element size, elx ~ 0.04 
mm and ely ~ 0.045, differed from the predictions obtained with the mesh illustrated in 
Figure 6 by less than 3%, and hence the solution was considered to have converged. 
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Figure 6. Mesh size and element selection 


BLIND PREDICTION RESULTS 


The sequence of events observed experimentally consists of an initial matrix-crack 
that develops between interfaces 1 and 2 within the taper region. That crack 
transitions into a delamination that grows along “interface 1” and migrates up via a 
matrix-crack to “interface 2”. This sequence of events was correctly simulated as 
illustrated in Figure 7. The main difference between simulated and observed damage 
events was the initial matrix-crack location, which was predicted near interface 2, at 
the weak singularity, in the taper region. The simulations also predicted multiple 


migration attempts prior to the final migration. Although these could not be observed 
at the edge of the specimens tested, X-ray images of tested specimens revealed 
migration attempts at the center of the specimens [1]. 

The blind-predictions for static loading are shown in Figure 8, together with the 
experimental load-displacement data. Two simulations, corresponding to assumed 
transverse tensile material strengths of Yr = 64 and 127 MPa, are depicted by red and 
black lines, respectively. These two values define lower and upper bounds of material 
strength values as measured by 90’ tensile or 3-point bend tests, respectively [13, 14]. 
As can be observed, the simulations showed a strong dependence on the assumed Yt. 
In both cases, damage initiation is predicted in the 90" ply near interface 2, at the weak 
singularity. 
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Figure 7. Final simulated damage state, obtained with Yr =127 MPa. Red represents elements 
completely opened. 
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Figure 8. Comparison of the load-displacement experimental response to simulations obtained with 
= 64 and Yr= 127 MPa. 


The blind-predictions for migration location are compared with experimental data 
in Figure 9. Results from blind predictions are shown for the two transverse strength 


values of Yr = 64 and 127 MPa. Migration location is measured from the load 
application point, as depicted in Figure 7. Simulations show that migration location 
occurs at greater distance from the load application point for Yr = 64 compared to the 
results obtained assuming Yr = 127 MPa. 
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Figure 9. Comparison of migration location obtained from the experiments and the simulated 
migration location obtained with Yr= 64 and Yr= 127 MPa. 


Overall, the results correlate well with the experiments. The peak load prediction 
obtained with higher Yr was closer to that experimentally obtained. The prediction 
obtained with Yr= 64 resulted in a significantly lower peak load prediction. However, 
in the taper region, the “as manufactured geometry” was different from the geometry 
idealized in the current analysis, as highlighted in Figure 10. Hence, the better 
agreement alone may not necessarily indicate Yr = 127 MPa is a more adequate 
transverse strength value. The migration process was also seen to be affected by the 
transverse tensile strength. In this case the migration distance corresponding to Yt = 
64 MPa was closer to the average migration distance measured experimentally. 
However, the migration distance obtained with Yr = 127 MPa was also within the 
experimental scatter. 
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Figure 10. Qualitative comparison between the “as manufactured” and the “as designed" ply location 
near the taper region. 


SENSITIVITY STUDY: EFFECT OF Yr 


In this section, the effect of strength on the peak load and migration location are 
further investigated and discussed. To do so, an additional simulation with Yr = 32 
MPa was performed, and the results obtained compared with the results obtained 
previously for Yr= {64, 127} MPa. 

Figure 11 compares the peak load predictions obtained assuming Yr = {32, 64, 
127} MPa. As can be seen, the peak load is very mildly affected by Yr below 64 MPa. 
This indicates that for Yr lower than 64 MPa the onset is controlled by the fracture 
toughness. 
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Figure 11. Variation in peak load prediction with Yr. 


The effect of Yr on migration location is particularly interesting, since crack 
deflection or kinking is often treated in brittle materials using Linear Elastic Fracture 


Mechanics (LEFM) based kinking criteria [18, 19]. In Figure 12, the migration 
location results obtained with Yr = {32, 64, 127} MPa are compared. The results 
show that despite the apparent similarity in peak load, Yr = 32 and Yr = 64 MPa lead 
to significant differences in migration location prediction. Additionally, while a lower 
Yt = 32 MPa leads to a conservative prediction in terms of peak load, it provides a 
non-conservative prediction in terms of migration location. The simulated damage 
state obtained with Yr = 32 MPa and Yr = 127 MPa, Figures 7 and 13, respectively, 
are also markedly different. With Yr = 32 MPa, prior to migration, a region of 
partially opened elements, is obtained, whereas with Yr = 127 MPa, only a few 
isolated migration attempts are observed. This finding indicates that the delay in crack 
migration is associated with the development of a migration attempt region with 
partially opened elements, akin to a process zone, prior to migration. The development 
of this process zone, composed of multiple partially opened elements, is originated by 
the reduction of strength while the fracture toughness remains constant. This is 
revealed in the present approach by not limiting the number of matrix-cracks and/or 
their spacing. In doing so, the methodology enables a more realistic material response 
and provides additional insight regarding the adequacy of the material 
properties/parameters used. 
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Figure 13. Final simulated damage state, obtained with Yr = 32 MPa. Red represents elements 
completely opened. 


SUMMARY 


As part of the NASA Advanced Composites Project (ACP), the clamped tapered 
beam specimen has been designed to provide validation data for progressive damage 
analysis models. The present work focused on the results obtained with a 3D discrete 
crack approach based on the Floating Node Method. The approach used provided 
adequate agreement with experimental measurements and shows potential in 
accurately simulating matrix-crack/delamination interaction. Three key components 
are identified in order to be able to model matrix-crack/delamination interaction: 

e accurately represent the kinematics of the matrix-crack/delamination 
interaction 

e model the matrix-crack angle and its through-thickness mode I 
propagation 

e recognize delamination and matrix-crack interaction are part of the same 
damage process and model them as such, rather than as separate 
independent events/damage modes. 

While fair qualitative agreement may be obtained with different degrees of 
approximation to any of the three components highlighted above, e.g. [4, 19, 20, 21, 
22], an accurate predictive strategy should aim at considering all of them. In the 
methodology used, of the three components identified above, the approximation used 
when modeling angled matrix-cracks and subsequent propagation through-thickness is 
recognized as the main aspect requiring further evaluation. 

The blind predictions also revealed challenges associated with obtaining 
characterization data for high-fidelity meso-scale modeling frameworks using 
cohesive zone crack approaches. Particular attention was given to YT as a result of its 
effect on peak-load predictions. Results suggest that the choice of Yr may not only 
affect the peak load predictions but also, and even more significantly, affect matrix- 
crack delamination interaction, which in turn may lead to an incorrect simulation of 
progressive damage. 

A wide range of values for Yr can be found in literature, as well as evidence of its 
volume dependency [16], which may undermine the experimental determination of a 
single value for Yr. In the present work, values of Yr leading to good qualitative and 
quantitative response were found to be high (~ 127 MPa) and pose challenges 
regarding the size of elements needed when using cohesive zone approaches, and 
hence its practical application. 

While lower values of Yr may provide conservative estimates of initial crack 
onset, they may also lead to incorrect simulation of crack interaction and hence 
erroneous post-peak predictions. This observation suggests that the commonly used 
approach of lowering the strength while keeping the fracture toughness constant, to 
enable the use of coarser meshes when simulating damage (in particular delamination) 
with cohesive zone elements, is not recommended in the predictive simulation of any 
component where crack initiation, small crack growth, and multiple crack interaction 
may occur. 
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